%write annual harmonics of ocean velocity
%as input to alexei's program
load d:\manoj\ecco216\daily_IV;
load d:\manoj\ecco216\daily_IU;
L = fday_ecco > datenum(1900,1,1);

mean_IU = squeeze(nanmean(IU(L,:,:),1));
mean_IV = squeeze(nanmean(IV(L,:,:),1));

load c:\manoj\ocean\ECCO_Cord;
[X1,Y1] = meshgrid(V_long,U_lat);

IUr = griddata(U_long,U_lat,real(mean_IU),X1,Y1,'linear');
IVr = griddata(V_long,V_lat,real(mean_IV),X1,Y1,'linear');

IUr(:,360) = IUr(:,1); %to remove NaN's from gridding
IVr(160,:) = 0;

IUr_n = [zeros([10,360]);IUr;zeros([10,360])];%padding zeros in north & south
IUi_n = zeros(size(IUr_n));
IVr_n = [zeros([10,360]);IVr;zeros([10,360])];
IVi_n = zeros(size(IVr_n));

data_iu = complex(IUr_n,IUi_n);
data_iv = complex(IVr_n,IVi_n);

load C:\Manoj\projects\tsunami\Bz;
Bz = rot90(Bz);
data_ur = flipud(real(data_iu));
data_ui = flipud(imag(data_iu));
data_vr = flipud(real(data_iv));
data_vi = flipud(imag(data_iv));
%all the date above are from 0.5 deg to 359.5 deg longitude


% prd = [4000]; % The mean flow has very high period -> 4000 days !
% 
% fp = fopen('d:\manoj\temp\a.project','wt');

  
CUr = -data_vr.*Bz.*1e-9.*3.2;
CUi = -data_vi.*Bz.*1e-9.*3.2;
CVr = -data_ur.*Bz.*1e-9.*3.2;
CVi = -data_ui.*Bz.*1e-9.*3.2;

%29..11..2007 I found that B_phi and B_theta (phi in particular) showing
%strange (Bp > 5nT at 400km) magntitude and pattern. Br is OK
% This also is present in the annual harmonics computations
%To check whether this is due to the calculation of impressed current
% I tried a series of combinations as shown below. but they made completely
%absurd results. also tried to "circshift" the bz and velocity data
%(idea was there may be an offset wrt the the conductance data in model
%file . that also failed.

%  CUr = -data_vr.*Bz.*1e-9.*3.2;%eccom_test WRONG WAY
%  CUi = -data_vi.*Bz.*1e-9.*3.2;% this produces absurd values
%  CVr = data_ur.*Bz.*1e-9.*3.2;% and patterns (eccom_test.res,
%  eccom_test.source)
%  CVi = data_ui.*Bz.*1e-9.*3.2;

%  CUr = data_vr.*Bz.*1e-9.*3.2;%eccom_test1
%  CUi = data_vi.*Bz.*1e-9.*3.2;%WRONG
%  CVr = -data_ur.*Bz.*1e-9.*3.2;%Produces wrong amplitude and patterns
%  CVi = -data_ui.*Bz.*1e-9.*3.2;
%  CUr = data_vr.*Bz.*1e-9.*3.2;%eccom_test2
%  CUi = data_vi.*Bz.*1e-9.*3.2;%WRONG
%  CVr = data_ur.*Bz.*1e-9.*3.2;%s
%  CVi = data_ui.*Bz.*1e-9.*3.2;
% 

%CU = CU*1000; % Compression for e format 20 oct 2005
%CV = CV*1000; % This multiplication is compensated by scale factor in source file
%fprintf(fp,'%s\n',['    ' sprintf('%10.5f',prd) '     ' sprintf('tsu%d',i) ' 1x1_mantle_k3.model ' sprintf('tsu%d',i) '.source' ' 1']);
fid = fopen('c:\manoj\ocean\temp\eccom.source','wt');
fid1 = fopen('c:\manoj\ocean\x3dg\header.txt','rt');
%the following format for alexei's x3dg is tested ok 20 oct 05
           for ii = 1:34,
                s=fgetl(fid1);
                fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
       fprintf(fid,'%12.4e\n',CUr(in,jn));
   end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
        fprintf(fid,'%12.4e\n',CUi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',CVr(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
  
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',CVi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:37,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
            fclose(fid);
            fclose(fid1);

fclose all;